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Summary 

This report describes a computer program for calcu- 
lating the flow field and flame propagation in a turbulent 
combustion tunnel. The program employs an algorithm 
for turbulent combustion modeling described by 
Ghoniem, Chorin, and Oppenheim (refs. 1 and 2). It uses 
the random vortex method of Chorin (ref. 3), which has 
the advantage of allowing the turbulent field to evolve as 
a fundamental solution of the Navier-Stokes equations 
without averaging or closure modeling. The flame front, 
which is treated as an interface between the reactants and 
the products, is computed by using an algorithm 
developed by Chorin (ref. 4) and exploiting a technique 
by Noh and Woodward (ref. 5) for simple line interface 
calculation. The flow is treated as “slightly 
compressible.” This allows for the increase in specific 
volume in the course of combustion to modify the flow 
field without any acoustic effects in order to model the 
effect of significant changes in density that occur at 
relatively low Mach numbers. 

The program provides numerical and graphical infor- 
mation about the vorticity field, the flame front location, 
and the composition profiles expressed in terms of the 
ratio of reactants to products, all as time-evolving fields. 
These data can be used to obtain average and stationary 
statistics of the turbulent field. In addition, overall 
performance parameters like the turbulent flame speed 
and the reattachment length behind a rearward-facing 
step in a channel can be deduced from these data. 

The program is rendered some flexibility by making 
most of its operation essentially independent of the 
geometry of the enclosure containing the flow field. 
Although most calculations are carried out in a universal 
transformed plane, the description of the geometry of the 
channel is restricted to a group of subroutines in order to 
allow a variety of configurations to be treated with a 
minimum of modifications. In this report the geometry 
of a channel with a rearward-facing step of an arbitrary 
height is used as an example. 


Introduction 

Turbulent combustion involves a complex set of inter- 
actions between essentially unstable fluid flow 
phenomena and transient exothermic chemical reactions 
that produce a significant increase in the specific volume 
of the flowing mixture. In general, these interactions 
influence the structure of the flow field as well as the 
progress of chemical processes and lead to a highly 
untractable situation, especially when basic mechanisms 
of both phenomena represent substantial challenges in 
their analyses. 


On one hand, chemical heat release, accompanied by 
large and rapid changes in specific volume, results in high 
flow acceleration that modifies and augments the tur- 
bulent field. It has been experimentally shown that when 
flame advection by the flow field is significantly larger 
than flame advancement due to normal burning speed, 
different rates of energy release can trigger different 
modes of instability and turbulent eddy interaction in the 
flow field according to Keller et al. (ref. 6). On the other 
hand, combustion fronts are essentially convected by the 
flow field, and their geometries are determined by the 
flow field structure. Thus overall reaction rates, 
controlled by mixing hot products and cold reactants, are 
strongly influenced by the turbulent field. The closed 
loop of interaction thus formed between turbulence and 
combustion makes it necessary for models of these 
processes to take their mutual influence clearly into 
account. Mellor (ref. 7) and Mellor and Ferguson (ref. 8) 
demonstrate that this complex interaction hinders 
numerical models from reproducing experimental data; 
Bray (refs. 9 and 10) points out the difficulties of 
accurately incorporating such interactions into 
computational schemes. 

In principle, the problem can be described in terms of a 
set of conservation equations coupled with a chemical 
kinetic scheme for a particular reaction as shown, for 
example, by Williams (ref. 11). In many practical 
situations one cannot expect a viable solution for this 
system of equations since they are highly nonlinear and 
intrinsically unstable and their fundamental solutions 
vary over a wide range of time scales. Moreover, in most 
cases of interest the equations are elliptic, with recircu- 
lating fields and sharp gradients in the solution domain, 
and thus contain a wide spectrum of length scales. 

Both the initial and boundary conditions are statistical, 
and the solutions are expected to be very sensitive to 
small perturbations in these conditions according to 
Libby and Williams (ref. 12). The effect of these 
perturbations cannot be played down. In fact, 
instabilities introduced in the flow field by the growth of 
perturbations can, paradoxically enough, enhance the 
stability of turbulent flames. This has been demonstrated 
by the numerical analysis developed in reference 2 to 
model the fluid mechanical phenomena of turbulent flow 
in a combustion tunnel. 

The model presented herein establishes the link 
between turbulence and combustion in terms of the 
primary dynamic variable: the velocity field. Turbulence 
is solved from first principles in terms of the dynamics of 
the vorticity field by using the random vortex method. 
The random vortex method computes the. vorticity field 
as an equivalent number of discrete elements, called 
vortex blobs. The interaction of the vortex blobs is 
determined, and they are moved with each time step. 
Diffusion is represented by the random walk of the blobs. 


This method neither imposes an averaging procedure, 
which in turn requires closure modeling, nor introduces 
extra numerical diffusion into the solution, which can 
dampen the inherent natural instabilities of the flow 
field. Thus it can reproduce detailed turbulent eddy 
interactions in rapidly changing fields. The mechanical 
effects of combustion are modeled by the motion of the 
flame front, which is composed of advection and normal 
propagation as a laminar flame, as well as by a 
distribution of volumetric sources inside the combustion 
zone that introduces the effect of combustion 
exothermicity into the velocity field. The linkage between 
turbulence and combustion in this model provides the 
mechanism through which the interaction between the 
two fields evolves with time and spreads throughout the 
field. Thus the algorithm is particularly useful for the 
study of such phenomena as inflammation, flame 
instabilities, and transient flow fields occurring in 
internal combustion engines. 

The report is divided into three sections: the model, the 
algorithm, and the program. Under the model the 
equations used in formulating the problem are presented 
along with an overall outline of the solution steps. Under 
the algorithm the solution procedures are described in 
detail and the formulas used to implement each step are 
given. Under the program the organization of various 
subroutines in the code, as well as the function of each, is 
presented. The computer program, which is available 
from COSMIC, University of Georgia, Athens, Ga. 
30602, contains an extensive list of variable definitions 
and comments on major steps to make it readily 
understandable. Examples of computer graphics motion 
pictures of the vorticity field and the flame front 
interface are presented in film supplement C-308, which 
is available on loan and described at the back of this 
report. 

The authors are grateful to Professor Alexandre J. 
Chorin for his help and advice throughout the course of 
this work. He also provided the original code of the flame 
propagation algorithm. Mr. Youwei Dai produced part 
of the computations presented in this report. 

Symbols 

d influence factor of vortex sheet 

F differential transformation function, d^/dz 

f fractional volume of burned medium in cell 
H\ height of incoming flow channel 

H 2 height of main flow channel 

h length of vortex sheet 

he side length of cell 


/ Vm 

Ns number of source blobs 

Ny number of vortex blobs 

n surface normeil vector 

R Reynolds number 

r position vector, x,y 

S burning speed 

normal burning speed 
t time 

A/ time step 

Ate combustion time step 

u velocity in x direction 

u velocity vector in physical plane, u,v 

V velocity in y direction 

w complex velocity, u — iv 

x,y Cartesian space coordinates 

z complex coordinate in physical plane, x+iy 

a expansion ratio, Fl^/Flx 

r circulation of vortex blob 

7 circulation per unit length of vortex sheet 
A source blob strength 

h thickness of sheet layer 

e volume source strength 

f complex vector in transformed plane 

77 Gaussian random number 

y dynamic viscosity of fresh mixture 

p specific volume ratio across flame 

^ vorticity 

p density of fresh mixture 

Po blob core radius 

Subscripts: 

b burned medium 

c combustion 

/ flame front 

ij for /th blob summed over all other j blobs 
s source velocity 

V vortex velocity 

6 produced by combustion 

b evaluated at edge of vortex sheet 

f produced by turbulence 

0 reference point 

00 incoming flow 

Superscripts: 

— complex conjugate 

* sheet property 
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The Model 

The analytical model is described schematically by the 
block diagram shown in figure 1. Each block contains 
either the relevant equations used to model a physical 
process or the link between two subsequent processes. 
The diagram is constructed of two loops that are linked 
together by the total velocity field u. The loop on the left 
side shows how vortex dynamics is applied to solve the 
turbulent field by the random vortex method (refs. 3, 13, 
and 14). The key element in this loop is the vorticity 
which is treated as a set of discrete elements that are 
transported by convection due to the total velocity field u 
and by diffusion expressed in terms of a random walk. 
Vorticity is created on solid boundaries to satisfy the no- 
slip condition and produces the solenoidal part of the 
velocity field u^. 

The loop on the right side of figure 1 illustrates how 
combustion is incorporated into the model. The flame 
front, at locations specified by the vector tp is treated 
locally as an interface between two incompressible media: 
the reactants and the products. The front is advected by 
the total velocity field u and propagates in the direction 
normal to its surface at a prescribed normal burning 
speed The location of the interface, as well as its 
displacement, is calculated by using an algorithm that 
utilizes the concentration field scalar /, which is defined 
as the local fractional volume of burned gas (ref. 4), to 
describe the composition field. The scalar field / is 
determined for a grid located in the flow. The changes in 
/ due to burning are employed to evaluate the specific 
volume source strength €, which is the key element of the 
combustion algorithm loop. These sources produce the 
irrotational part of the velocity field u. 


A typical graphical output that shows various 
computational elements is presented in figure 2. The dots 
on the tails of the lines represent the positions of the 
vortex blobs, and the length and direction of the lines 
represent the magnitude and direction of the velocity. 
The dashes represent the edge of the flame. 

The algorithm is unsteady and forward marching in 
time. The solution procedure is as follows: 

(1) Start with potential flow 

(2) Add turbulence at wall (vortex sheets) 

(3) Do random walk 

(4) Find flame 

(5) Move with flow 

(6) Move with laminar flame speed 

(7) Repeat, starting with step (2) 

Steps (1) to (3) are flow calculations; steps (4) to (6) are 
flame calculations. The solution is begun at / = 0. There is 
potential flow throughout the channel except on the 
boundaries where vortex sheets are added to satisfy the 
no-slip condition. The flow calculations are performed 
first and then the flame calculations. Computations 



^Vortex sheets in 
Iwundary layer 


Figure 2. - Graphical output of MIMOC. 


Vortex dynamics 


Flame propagation 


Solenoidal field 


Kinematics 


Irrotational field 



Boundary diffusion 


Burning 


Figure 1. - Structure of MIMOC algorithm. 
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continue until enough iterations are reached to obtain 
steady statistics. 

To calculate the vortex dynamics, the flow field is 
partitioned into a core, which represents the interior 
region of the flow field, where vortex blobs are located, 
and a numerical boundary layer, or several of them if 
more than one wall is considered, where vortex sheets are 
located. The different operations of the random vortex 
method and the relationship between the two kinds of 
vortex elements, namely vortex sheets and vortex blobs, 
are described next. 

Each boundary layer extends over a very small 
d istance 6, normal to the wall, in the order of magnitude 
V2 Ar/iR, where At is the time step and R is the Reynolds 
number. Within a boundary layer, vortex elements, gen- 
erated on the wall to satisfy the no-slip boundary 
condition, are discretized along the wall direction into 
vortex sheets with a finite length h. The influence of these 
vortex sheets is restricted locally to a domain normal to 
their direction, the wall direction, and their diffusion is 
taken only normal to the wall. These elements provide a 
fine resolution in the near-wall region, and they satisfy 
the conservation of vorticity by virtue of their reflective 
properties (ref. 13). Their use also significantly reduces 
computational time as a consequence of their confined 
local influence in regions of high vorticity concentration. 

In the interior, vortex elements become vortex blobs as 
they diffuse out from numerical boundary layers. The 
induced velocity field of a vortex blob is fully symmetric 
although the blob movement by diffusion is completely 
unbiased. The transition between sheets and blobs is 
smoothly maintained in the solution of the vorticity 
transport equation by providing a region of overlap of 
the two elements to allow one element to transform to 
another while it maintains the value of its circulation and 
identity. The continuity of the velocity field is maintained 
by extending the velocity induced by vortex blobs into the 
vortex sheet region. 

Because the combustion model is based on the 
assumption of fast chemistry similar to Williams (ref. 
15), either the reactants or the products may exist at the 
same place with an infinitely thin interface separating 
them. The conversion of one component of the mixture 
into the other is controlled by the combined effect of two 
factors: (1) the normal burning speed which accounts 
for molecular diffusion processes as well as chemical 
reaction kinetics, and (2) the kinematics of the interface, 
governed by the distribution of the concentration field 
/(r). Besides, the interface is transported, or advected, by 
the flow field while it modifies that field by the velocity 
produced because of the volumetric expansion across the 
flame. 

The concentration field is used as a way of digitizing 
the geometry of the flame. The plane is divided into a 
number of square cells by using a grid of a fixed mesh 


size, and / is given a value in each of these cells to 
describe the mixture distribution. The interface ifiside a 
cell is treated as a simple straight line whose location is 
determined by the value of /. Its orientation as either a 
horizontal or a vertical line depends on the value of /in 
the cell and its four neighboring cells. This technique is 
described in more detail in the next section. The motion 
of the interface line is determined by its speed, and the 
concomitant displacement of the fluid on both sides is 
used to change the value of /in the cells on both sides. 
When the direction of the speed is not known, like the 
propagation of the flame normal to itself, the scheme 
suggested by Chorin (ref. 4) is applied. According to 
Chorin, the advancement of the flame due to burning can 
be evaluated by Huygen’s principle on recognizing that 
its front, propagating in the direction normal to its 
surface, forms an envelope to its displacement in all 
possible directions. 

By applying the continuity requirements reference 2 
used this algorithm to calculate the amount of reactants 
converted into products at each time step and 
consequently the volumetric change of the flow field due 
to combustion. This change gives rise to an irrotational 
velocity field, according to Batchelor (ref. 16) and Chorin 
and Marsden (ref. 17), which is added to the vortex 
velocity field to furnish a complete solution of the 
turbulent combustion system under investigation. 


The Algorithm 

A pragmatic description of the algorithm is presented 
in this section. It includes all of the formulas used for the 
computations but does not provide their derivation. The 
details of the theoretical background can be found in 
reference 2. All quantities used in the program are 
nondimensionalized with respect to the reference param- 
eters. Length is normalized with respect to the height of 
the main channel H 2 . Velocities are normalized with 
respect to the incoming flow velocity Uc», and time is 
normalized with respect to //qo/Uoo- Thus the Reynolds 
number is defined as PU 00 // 2 /M, where p is the fresh 
mixture density and p is its dynamic viscosity. 

The computational algorithm is divided into two 
sections: computations of vortex dynamics and compu- 
tations of combustion. In the first section, vortex 
elements in the interior of the flow, called vortex blobs, 
are moved in both the physical z-plane and the trans- 
formed upper-half f plane. The transformation is 
accomplished by the general Schwarz-Christoffel trans- 
formation, which allows the solution to be mapped into a 
variety of physical geometries. This technique results in a 
computer program with most of the subroutines being 
geometry independent. The transformed plane is used to 
facilitate the computations of a velocity field with a zero 


4 



normal component along the boundaries. Vortex blobs 
are propelled with their own induced velocity as well as 
with the velocity field imposed by both the incoming 
potential flow and the volumetric source blobs. 

The physical plane and the transformed plane for a 
channel with a rearward-facing step are shown in figure 3. 
The positions of vortex elements on the boundary walls, 
called vortex sheets, are then updated after calculating 
the velocities of the sheets from the velocity field imposed 
on them by the interior field. If the no-slip condition on 
the walls is not satisfied, new sheets are created and 
added to the existing vorticity field. In the second section 
of the algorithm, combustion calculations are performed 
by advancing the flame in two fractional computational 
steps. In the first fractional step the flame is advected by 
the velocity field; in the second step it is displaced by the 
normal burning speed of the corresponding laminar 
flame in order to account for combustion. The strength 
of volumetric sources is computed according to the 
combined effect of the burning speed and the specific 
volume ratio across the flame. Their location is 
established by the position of the flame front. 

Computations of Vortex Dynamics 

The computations start at time r = 0, with the flow field 
produced by the incoming potential flow only. The 
tangential velocity produced by this flow along the solid 
walls is computed and used to produce the first vortex 
sheet in the field, with the proper strength to cancel this 
velocity. The vortex sheet is divided into a number of 
sheets of length h along each wall, and each of them is 
subdivided into a stack of smaller sheets normal to that 
wall. The diffusion of these stacks away from the walls 
and into the interior field is the source of vorticity in the 
flow during the course of the computations. 



Hi] Ueo:^ 
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(a) Physical plane. 

(b) Transformed plane. 

Figure 3. - Geometry of physical and transformed planes. 


The computation for each new time step, ^ + Af, starts 
by updating the vorticity field Hx,y). This is accomp- 
lished in two steps: (1) translation of the vorticity field, 
and (2) augmentation of the vorticity field. The 
computation of each step is presented in detail. 

Translation of vorticity /fe/d. -New coordinates of 
each vortex blob are calculated from 

z{t+At) =z(0 +» (1) 

where z =x+iy is the complex coordinate of the center of 
a vortex blob in the z plane, u = « -h /y is the velocity in the 
complex plane, and v-Vx~^^Vy is the two-dimensional 
random walk modeling the diffusion of vorticity of which 
each component is computed as a Gaussian random 
number with zero mean and variance 2 At/R. (For more 
details about modeling diffusion by random walk, see 
Ghoniem and Oppenheim, ref. 18.) 

The total velocity in the z plane is evaluated by first 
calculating the velocity in the f plane as a result of the 
following effects: 

(1) Incoming flow w,- 

(2) Vortex blobs and their images, the latter to satisfy 
the zero normal velocity condition along the solid walls, 

(3) Source blobs and their images, the latter to cancel 
their induced normal velocity along solid boundaries, 

The value of the velocity, expressed in terms of a complex 
velocity, can be written as 

Ny Ns 

M'(f)=W,(0+ ^ D W.iUj) (2) 

where 

H’,(D=“ (3) 


The velocity field induced by each individual blob, 
evaluated as the fundamental solution of the solenoidal 
field and irrotational field equations shown on the top 
corners of figure 1, is given by 

wAUj) = -'T/2^ Max[|r- f>oKf- fy) 

+ iTj ^ (4) 

•'2irMax[|r-fy|,Po](f-fy) 

for a vortex blob and by 
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Max[|r- fy) 

^ 

" 27 t Max[|r-fy|,poKf-fy) 


for a source blob. In these equations, f is the coordinate 
in the transformed f plane, f=conj(f), Vj is the circu- 
lation of a vortex blob, Ny is the total number of vortex 
blobs, Aj is the strength of a volumetric source blob, A^5is 
the total number of source blobs, and pq is the radius of 
the blob core. The absolute value of the vortex-induced 
velocity inside this core is taken as a constant. The radius 
of the core in the f plane is calculated by using the 
transformation 

Po(fy) = Po^^(fy)| (6) 

where pQ^h/ir is the radius of the blob core in the z plane 
and h is the length of the vortex sheet. The function F(f), 
the differential transformation that maps the z plane into 
the f plane, is defined as 


(7) 


where z is nondimensionalized with respect to the channel 
height. The function F used for the rearward-facing step 
with an expansion ratio a is 

/r-a2\l/2 

The velocity in the z plane is then obtained by the inverse 
transformation 

>v(z)=iva)F(f) (8) 

and u = conj[w(z)]. 

Because some vortex blobs may move into the 
boundary layers, the vorticity between the interior and 
the various boundaries must be resorted to make sure 
that the contribution of a vortex element is correctly 
taken into account in the next steps. The location of 
vortex blobs in the f plane is required in the next time step 
in order to calculate their induced velocity according to 
the above procedure. The location is obtained by 
applying the following formula: 

f(r-hAO = f(0 + t F(f)rfz (9) 

J Az 


where Az=z(/+A/) -z{t). The integration from the z 
plane to the f plane is performed by using a fourth-order 
Runge-Kutta method. 

Augmentation of vorticity field. -PiAAiiiondl vortex 
blobs appear in the interior field as a result of vortex 
sheets moving outside their numerical boundary layers, 
as shown in figure 4. These sheets are originally generated 
on the wall to cancel the tangential velocity on it. Outside 
a layer in the order of magnitude of the variance of the 
random walk, vortex sheets become vortex blobs with the 
same circulation in order to satisfy the conservation of 
the total circulation. Since vortex sheets can be reflected 
back into the field if they move outside it, as shown by 
Chorin (ref. 13), this layer helps to reduce the possibility 
of losing vortex blobs. 

To compute the position of vortex sheets, to transform 
sheets into blobs when necessary, and to generate new 
sheets on the walls, the following steps are taken for each 
wall: 

(1) The velocity produced by the interior field is 
computed on the walls. 

(2) The velocity is transformed into the wall 
coordinate system in which the wall is taken as the x axis 
and the normal to it as the 7 axis. This transformation is 
necessary to keep the calculations of vortex sheets 
independent of the wall orientation. 

(3) Vortex sheets on the wall are moved according to 
the formula 

Z/(/ + A/) = z/(/) +W/ At^^riy (10) 

where w/= w/H-Zu/ and 
5 

«/•= - 0.5 7, - E?/ dj (1 1) 

yi 

Here u^. is the velocity, calculated at x/, at the edge of the 
sheet layer produced by the interior field, 7/ is the circu- 
lation per unit length of the sheet, and dj is the influence 
factor defined as 



Figure 4. - Transformation of vortex elements near a solid wall between 
sheets and blobs. 
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W-=1- 
"y /I 

The summation is taken over all of the sheets that satisfy 
the following conditions: 

yi>yj 

\>dj>0 

Figure 5 provides a graphical interpretation of 
equation (11), where the zone of influence of a vortex 
sheet is marked by its shadow. The normal-to-the wall 
velocity component is 


and 

]pCi±{h/2)-Xi \ 

J h 

The summation is taken over all of the sheets for which 
and 

^/=MinO/,>>y) 

The y coordinate of the new location is checked for 
negative values; if >'/<0, the sheet is reflected back with 

yi= -yi. 

(4) The tangential velocity on the wall is computed by 
using a special case of equation (11). If y,=0, 



(12) 


Uv,,=y&{Xi) 


E T/rfy 

y = 0 


(13) 


where 




A Zone of dependence over point i 

B Zone of influence under sheet j 

C Zone of dependence around point i + L 

D Zone of dependence around point i - ^ 



If |wvv,^|^rniin> where r^in is an assigned value for the 
minimum circulation of a new vortex element, the 
appropriate number of vortex sheets are generated or 
removed on the wall. The numerical choice of these 
parameters is discussed in a later section. 

(5) Finally, the coordinates of all of the vortex sheets 
are checked for y>6, where 6 is the thickness of the 
numerical boundary layer. Vortex sheets that leave this 
layer are converted into vortex blobs in the interior field. 
Their coordinates are transformed back into physical 
coordinates and their locations in the f plane are 
calculated by using the appropriate modification of 
equation (9) 

f=fo+r (14) 


where Zo represents a reference point with a predefined 
transformation fo- The analytical transformation that is 
used to evaluate zo from fo is 

(log|^--log^i^) (15) 

where 

/cx2_f\l/2 

( w) 

Computations of Combustion 

The computations of combustion consist of two 
components: (1) flame advancement, and (2) dynamic 
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effects of combustion due to volumetric expansion. The 
flame motion is composed of advection and combustion. 
They are taken into account by fractional steps as 
follows. 

Advection. - A grid is first constructed in the physical 
plane to keep track of the volumetric concentration of 
products f in the combustion zone. The grid divides the 
plane into a number of cells of size as shown in fig- 
ure 2. So that the advection velocity can be calculated on 
the grid, the centers of these cells are transformed onto 
the f plane by using equation (14). Equations (2) to (5) 
are applied in each step to evaluate the velocity at these 
centers, which is then linearly interpolated to compute 
the velocity component required on the sides of the cells, 
as shown in figure 6(a). 

The orientation of the interface inside a cell is 
determined by comparing the value of /in the cell and its 
value at the four neighboring cells. If /=0 or /= 1, the 
cell is filled with reactants or products, respectively, and 
no computations are required in this cell. If 0</<l, 
the interface is constructed as one of the following 
possibilities: 

(1) Vertical neck: 

^•-u=0 

and 

X+ij=o 

(2) Horizontal interface: 
and 


and 

fij+\ = ^ 

Figure 6(b) shows the geometrical interpretation of these 
arrangements. To simplify the computations, motion in 
two directions is split into two fractional steps, each in 
one dimension. 

After the interface is located, its velocity is evaluated 
by interpolating between the velocities on both sides of 
the cell. The interface velocity is used to move the 
interface. The fluid on both of its sides is displaced in or 
out of the surrounding cells to accommodate this motion. 

Combustion. -TYic motion of the flame in the 
direction normal to itself due to combustion is calculated 
in the exact manner as advection but by using the normal 



(3) Vertical interface: 
and 

(4) Corner: 



II 


Vertical neck 

(a) Velocity components. 

(b) Interface geometries. 

Figure 6. - Computational elements for flame advection. 



Horizontal interface 
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burning speed 5^ in all possible directions of propagation 
instead of the advection velocity. Since the direction of 
propagation of the flame due to combustion is not known 
apriori, is given all possible directions and the actual 
flame motion corresponding to all of the directions is 
calculated. Thus the flame speed is taken as 


combustion step A/^ are used to calculate the strength of 
the volumetric sources according to 


A = 0.5(;--l)/j2^ 


(19) 


Sf= (Su cos $i, Sy sin Si) (17) 

where 


where p is the specific volume ratio across the flame. For 
each cell with A/c>A/nun, where A/^in is a negligible 
change in volume, a source blob is generated and placed 
at the center of this cell. The core size of this blob is 


' 4 



( 20 ) 


and 

/= 1 , 2 , . . ., 8 


where Us=0.5iv- 1)5„ is the maximum velocity generated 
by the flame propagation. 


According to Huygen’s principle, the actual flame 
propagation is the maximum advancement in all 
directions. Thus / can be updated as 

/(/ + A/) = Max/' (18) 

0</<8 


where 

/®-/(0 

Combustion exothermicity. - Figure 7 explains the 
dynamic effect of burning due to the increase in specific 
volume in which a particle moves away from the flame 
front a distance of 2Us after burning. To account for this 
effect in the two-dimensional calculations of the velocity 
field, changes in f due to flame advancement in the 

tA 



Choice of Numerical Parameters 

In choosing the input parameters to the computations, 
a compromise is made between the accuracy of the results 
and the efficiency of the computations in terms of time 
and storage. 

To limit the number of vortex elements produced and 
subsequently used in the computations, the following 
values are recommended: 


= 0.2 

(21) 

min 

(22) 


This value of Fmin results in the generation of four sheets 
in the first time step. Thus it provides a reasonable 
accuracy in calculating of the diffusion of these sheets 
into the field. Since the error of the method is expected to 
be O(A0 (ref. 3), the time step for the nonreacting flow is 
taken as 

At^OA (23) 

In the combustion calculations a cell size 

hc = 0.05 (24) 

was found appropriate to achieve a fine resolution of the 
flame front. Note that, since the flame surface is located 
within the cell, the resolution of the flame surface is not 
limited by The combustion time step must satisfy the 
Courant condition 


At< 


he 

2(Umax + S„) 


(25) 
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to keep the calculations stable. In most cases, Atc<Aty 
and several combustion steps are performed for each 
vortex step by choosing At—K Ate, where K is an integer. 


The Program 

Figure 8 depicts the overall flow chart of the MIMOC 
program. For the sake of simplicity, only the major steps 
and subroutines are shown. The program listing contains 
an extensive array of comments and variable definitions 
to further explain it in detail. There are no input data for 
the program. All of the necessary parameters of the flow 
field are initialized in MIMOC; those for the walls are 
calculated in WLPNTS. If combustion computations are 
to be performed, subroutine INTCMB is called to initi- 
alize the combustion parameters such as the burning 
speed and the density ratio across the flame. If these 
parameters are changed, the program must be 
recompiled. The thick solid lines indicate that the control 
is in MIMOC. The thin solid lines identify the control 
with either of the two subroutines BOUNDR or 
COMBUS. The double box signifies the role subroutine 
POTVEL plays in the center of the calculations in the 
code. 

At each time step, t = t-\-At, subroutine BOUNDR is 
called four times, one for each wall in the channel. 
BOUNDR acts like a master program for the application 
of the no-slip boundary condition along solid walls. To 
accomplish that, it starts by calling POTVEL to calculate 


the velocity on the wall due to the interior field. 
Subroutine SHEETS is called to calculate the velocity of 
vortex sheets and the velocity at the wall by using 
subroutine SHEVEL. If the net velocity at the wall is 
larger than a minimum, vortex sheets are generated or 
removed. SHEETS is also used to move the existing 
vortex sheets by using the velocity produced by the 
interior field as well as their induced velocity. Finally, 
BOUNDR calls SHEBLB to check the position of the 
vortex sheets and to transform those that leave the sheet 
layer into vortex blobs. Subroutines NRMWL2 and 
PRLWL2 are used to transform the coordinates of vortex 
sheets from wall coordinates to physical coordinates 
when they become vortex blobs. These blobs are fed 
immediately to POTVEL to account properly for their 
induced field, as indicated by the short -dashed line. 

The control goes back to MIMOC after the 
calculations at the boundaries are finished. The position 
of the vortex blobs is updated by moving them according 
to their induced velocity and by random walk. MIMOC 
acts also as a sorter for vortex blobs in order to assign 
those near boundaries to the proper sheet layer. Vortex 
blobs that leave the interior zone are assigned back to 
SHEETS, as indicated by the long-dashed line, after their 
coordinates are adjusted in either NRMWLl OR 
PRLWLl. 

Subroutine COMBUS acts as a master program for the 
computations of flame propagation. First FLAME is 
called to move the flame by combustion. Here the normal 
burning speed S„ is used, and subroutine XNOH is called 
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Figure 8. - Flow chart of the MIMOC main program. 










to locate and move the interface inside each cell. This 
motion is employed in COMBUS to calculate the strength 
of the volumetric sources. COMBUS calls POTVEL to 
calculate the advection velocity at the center of the cells 
and uses it to calculate the velocity on the sides of these 
cells. This velocity is used in subroutine ADVECT to 
move the flame by calling XNOH for each cell. In 
MIMOC the results are written on a permanent file for 
later use and the time is updated again. 


Results 

Samples of the results obtained by using the program 
are shown in figures 9 to 12. The samples involve the 
vorticity field, the average velocity field, the turbulent 
statistics, the flame front geometry, and the 
concentration contours in a combustion channel with a 
rearward-facing step. (Similar results, computed for a 
free jet by using a modified version of the program, are 
presented in ref. 19.) The graphical outputs of the 
vorticity field and the flame front interface were plotted 
by employing computer graphics on a 16-mm black-and- 
white motion picture film. The concentration contours 
were displayed on a color image processor by assigning a 
relationship between color and the value of / and were 


recorded on color motion picture film. Examples of tht 
output are provided in the motion picture supplement 
(C-308) of this report. 

The computational time required to run this code can 
be long, as much as 4 to 5 hours on an IBM 370. As the 
number of vortex blobs increases, the running time per 
computational step grows as O(N^). However, in con- 
fined flows, such as the one considered in this report, Ny 
remains approximately constant. The number of vortex 
blobs used depends on the length of the vortex sheets h, 
the time step A/, and the circulation of the vortex blobs 
Tjnin- The number of blobs determines the accuracy of 
the solution and the smallest scale of turbulence that can 
be resolved. 

All of the computations were performed for a 
Reynolds number of 10 000. Presented in figure 9 is the 
vorticity field in flow over a rearward-facing step with an 
expansion ratio of 2. Figure 9(a) shows the development 
of the flow field by presenting vortex velocity vector 
fields that trace the motion of all of the vortex blobs 
included in the solution at successive dimensionless time 
intervals, each equal to 50 computational steps of 0.1. 
Figure 9(b) shows the growth of a large-scale eddy traced 
at time intervals equal to five computational steps. The 
formation of large-scale turbulent eddies, by the 
interaction of the incoming flow and the recirculating 
flow, can be readily seen in these results. 


Dimensionless 

time, 

t 




{a) Development of flow field. 

(b) Growth of large-scale eddy. 

Figure 9. - Vorticity field in channel with rearward-facing step. Expansion ratio, 2. 
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Figure 10 shows the average velocity profiles and the 
turbulent stress profiles in the same channel but with a 
smaller step (expansion ratio of 1.5). There are two sets 
of profiles for each variable. The set on the top cor- 
responds to averaging between time steps 291 to 348, the 
set on the bottom, between time steps 158 to 200. The two 
sets resemble each other qualitatively although the top set 
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(a) Axial velocity. 

(b) Cross -stream velocity.^ 

(c) Turbulent normal stress j^. 

(d) Turbulent shear stress -u'V. 

Figure 10. - Velocity profiles and turbulent statistics profiles. Reynolds 
number. 10 000; expansion ratio. 1.5. 


is closer to steady-state conditions than the bottom set. 
This comparison shows that at least 300 steps are needed 
to build valid statistics. According to the average 
streamwise velocity profiles, shown in figure 10(a), the 
channel can be divided into three regions. The first region 
is dominated by the recirculation bubble with strong 
negative velocities. The second region starts at the re- 
attachment zone with a boundary layer developing under 
a favorable pressure gradient. In the last region the two 
boundary layers, growing on both walls, converge to 
form the start of a fully developed flow in the channel. 
The average cross-stream velocity profiles are depicted in 
figure 10(b). Here the recirculation zone is marked by a 
positive value of this velocity component, but the rest of 
the channel is dominated by a negative vertical velocity 
that is due to the expansion of the flow to fill out the 
channel. The turbulent normal and turbulent shear 
stresses shown in figures 10(c) and (d), respectively, 
exhibit a maximum at the high shear zones, namely at the 
shear layer downstream from the step and in the 
boundary layers on both walls. 

Figure 11 displays the variation in the vorticity field 
and the flame front contour, shown as a demarcation line 
between cells with /=0 and those with />0, calculated 
for a channel with an expansion ratio of 2, = 0.02, and 

p = 4. These data correspond to a propane-air mixture 
burning at an equivalence ratio of 0.5. The figure consists 
of two sequential series of computer outputs. The 
sequence in figure 1 1(a) displays the process of ignition in 
the turbulent flow of figure 9 in a cell located on the 
centerline of the channel at a distance from the step equal 
to the height of the step. The sequence in figure 11(b) 
depicts the “stationary state” operation reached at 
/ = 26. 1 after ignition at the bottom corner at the moment 
when the medium was set in motion. The number of 
computational time steps, A/ = 0.05 each, between the 
solutions displayed here is 40 for figure 11(a) and 4 for 
figure 11(b). From these figures it can be concluded that 
the flame stabilizes itself on the outer edge of the large- 
scale eddies of the turbulent field and establishes an 
acceleration that reduces the length of the recirculation 
zone from that of a flow without combustion. These 
observations agree with the results of experiments by 
Ganji and Sawyer (ref. 20) and Pitz and Daily (ref. 21) in 
a similar channel with equivalent flow condition. 

A black-and-white image of a color photograph of the 
flow described above is shown in figure 12. The 
numerical solution is displayed in figure 12(a) and the 
experimental record, in figure 12(b). The color of the 
numerical results was assigned according to the value of 
/, producing areas of constant concentration of products; 
the color of the experimental results changed with the 
gradient of the local temperature. In the black-and-white 
image, the bright areas, in both cases, represent the zones 
of maximum burning. The agreement between them is 


Dimensionless 

lime. 
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(a) Ignition in turbulent flow of figure 9 In cell located on centerline of channel at distance from step equal to step height, 
(b) Stationary-state operation at t = 26. 1 after Ignition at bottom corner. 

Figure 11. - Flame geometry and vorticity field with rearward -facing step. Reynolds number. 10 000; expansion ratio. 2. 



(a) Numerical. lb) Experimental. 

Figure 12. - Concentration contours. 


remarkable, considering the simplifications used to 
produce the numerical results. A 16-mm motion picture 
supplement (C-308) shows the dynamic nature of the 
calculation. 


Conclusions 

A computer code for calculating unsteady, two- 
dimensional turbulent reacting flow is presented in this 
report. The model used for the analysis is based on the 
random vortex method of computing the turbulent flow 
field. Because the method is grid free and devoid of 
numerical diffusion, it is suitable for modeling the 
intrinsic physical properties of flows at large Reynolds 
numbers. It can analyze large gradients produced by the 
small eddies of turbulent motion, and it allows pertur- 
bations to grow into possible instabilities without 
damping or undue dissipation. The combustion of 
premixed gases is idealized by the propagation of a flame 
interface, located between reactants and products, due to 
convection and self-propagation. The expansion of the 
flow across the flame due to the exothermicity of the 
reaction is taken into account by using a distribution of 
volumetric sources within the burning zone. 

The program was used to study the flow field in a 
model combustor, formed by a rearward-facing step in a 
channel, in terms of the vorticity field, the velocity field, 
the turbulent shear stresses, the flame contours, and the 
concentration field. Results for the vorticity field reveal 
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the formation of large-scale eddy structures in the 
turbulent flow downstream from the step. Not only do 
these structures persist in the reacting case, but also the 
flame front almost surrounds the vorticity field and thus 
produces a situation in which the flame propagates along 
a streamline instead of normal to it. The process of 
stabilizing a turbulent flame behind a bluff body is 
illustrated by igniting a fully developed turbulent flow. 
The formation of a recirculation zone behind the step 
provides a continuous source of ignition by enhancing the 
contact between the products and the reactants. In 


addition, combustion is shown to reduce the length of the 
recirculation zone. The concentration field contours 
indicate that most burning occurs around the outer edges 
of the large eddies of the shear layer and that the 
combustion process extends far downstream in the 
channel. 


National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio, March 28, 1983 
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Appendix - Routines in Computer Program 


The 27 routines in the order that they appear in the 
computer program are listed here. 


Vorticity Field Routines 


MIMOC 

BOUNDR 

CFSCZ 

CRK4 

INIPUT 

NMTRM 

NRMWLl 


NRMWL2 


PLTFLM 

PNTWAL 

POTVEL 

PRLWLl 


main routine, initializes variables and 
directs program 

controls program flow to satisfy no-slip 
conditions at the walls 
differential transformation function 
Runge-Kutta integration 
starts integration procedure using ZXACT 
numerical transformation using reference 
point 

transforms from physical coordinates to 
wall coordinates for vertical walls when a 
blob becomes a sheet 

transforms from wall coordinates to 
physical coordinates for vertical walls 
when a sheet becomes a blob 
plots graphs 

checks on wall point transformation 
computes velocity in physical plane from 
points in complex plane 
same as NRMWLl except for horizontal 
walls 


PRLWL2 

RANDNM 

SHEBLB 

SHEETS 

SHEVEL 

TRNSPP 

WLPNTS 

ZXACT 


same as NRMWL2 except for horizontal 
walls 

generates Gaussian distribution of random 
numbers 

transforms vortex sheets into blobs 
adds or removes sheets to satisfy no-slip 
boundary condition 
calculates sheet velocity 
transforms point by using a reference point 
sets values of wall points where boundary 
conditions will be checked 
performs exact transformation from 
complex plane to physical plane 


Combustion Routines 


ADVECT 

COMBUS 

FLAME 

INTCMB 

XNOH 

WD 

WDP 


advances flame by advection 
controls program flow for combustion 
advances flame by laminar flame speed 
initializes combustion parameters 
moves individual flame interface 
generates uniform random numbers used 
by INTCMB and COMBUS 
same as WD but completely independent 


15 


References 


1. Ghoniem, A. F.; Chorin, A. J.; and Oppenheim, A. K.: Numerical 

Modeling of Turbulent Combustion in Premixed Gases. 18th 
Symposium (International) on Combustion, The Combustion 
Institute, 1981, pp. 1375-1383. 

2. Ghoniem, A. F.; Chorin, A. J.; and Oppenheim, A. K.: Numerical 

Modeling of Turbulent Flow in a Combustion Tunnel. Phil. Tran. 
R. Soc. London Ser. A, vol. 304, 1982, pp. 303-325. 

3. Chorin, A. J.: Numerical Study of Slightly Viscous Flow. J. Fluid 

Mech., vol. 57, pt. 4, 1973, pp. 785-796. 

4. Chorin, A. J.: Flame Advection and Propagation Algorithm. 

J. Comp. Phys., vol. 35, 1980, pp. 1-11. 

5. Noh, W. T.; and Woodward, P.: SLIC, Simple Line Interface 

Calculations. International Conference on Numerical Methods in 
Fluid Dynamics, 5th, A. I. Vooren and P. J. Zandbergen, eds.. 
Springer- Ver lag, 1976, pp. 330-339. 

6. Keller, J. O.; et al.; Mechanism of Instabilities in Turbulent 

Combustion Leading to Flash Back. AIAA J., vol. 20, no. 2, 1982, 
pp. 254-262. 

7. Mellor, A. M.: Turbulent-Combustion Interaction Models for 

Practical High Intensity Combustor. 17th Symposium 
(International) on Combustion. The Combustion Institute, 1978, 
pp. 377-387. 

8. Mellor, A. M.; and Ferguson, C. R.: Practical Problems in 

Turbulent Reacting Flows. Turbulent Reacting Flows, P. A. Libby 
and F. A. Williams, eds., Springer-Verlag, 1980, pp. 45-64. 

9. Bray, K. N. C.: The Interaction Between Turbulence and 

Combustion. 17th Symposium (International) on Combustion, 
The Combustion Institute, 1978, pp. 223-233. 

10. Bray, K. N. C.: Turbulent Flow with Premixed Reactants. 

Turbulent Reacting Flows, P. A. Libby and F. A. Williams, eds., 
Springer-Verlag, 1980, pp. 115-184. 

11. Williams, F. A.: Combustion Theory. Addison-Wesley, 1965. 


12. Libby, P. A.; and Williams, F. A.; Fundamental Aspects— 

Turbulent Reacting Flows. P. A. Libby and F. A. Williams, eds., 
Springer-Verlag, 1980, pp. 1-43. 

13. Chorin, A. J.: Vortex Sheet Approximation of Boundary Layers. 

J. Comp. Phys., vol. 27, 1978, pp. 428-442. 

14. Chorin, A. J.: Vortex Models and Boundary Layer Instabilities. 

SIAM J. Scientific Stat. Comp., vol. 1, no. 1, Mar. 1980, pp. 
1 - 21 , 

15. Williams, F. A.: A Review of Some Theoretical Consideration of 

Turbulent Flame Structure, Analytical and Numerical Methods 
for Investigation of Flow Fields with Chemical Reactions, 
Especially Related to Combustion, M. Barrer, ed., AGARD CP 
164, AGARD, 1975, pp. III-l to III-25. 

16. Batchelor, G. K.: An Introduction to Fluid Mechanics. Cambridge 

University Press, London, 1967. 

17. Chorin, A. J.; and Marsden, J. E.: A Mathematical Introduction to 

Fluid Mechanics. Springer-Verlag, 1979. 

18. Ghoniem, A. F.; and Oppenheim, A. K.: Random Element Method 

for the Numerical Modeling of Diffusional Processes. Presented at 
the 8th International Conference on Numerical Methods in Fluid 
Dynamics (Aachen, West Germany), June 28-July 2, 1982. 

19. Ghoniem, A. F.; Chen, D. Y.; and Oppenheim, A. K.: Formation 

and Inflammation of a Turbulent Jet. Presented at the Spring 
Technical Meeting of the Combustion Institute, Canadian Sec. 
(Banff, Alberta, Canada), May 9, 1982. 

20. Ganji, A. R.; and Sawyer, R. F.: An Experimental Study of the 

Flowfield of a Two-Dimensional Premixed Turbulent Flame. 
AIAA J., vol. 18, no. 7, July 1980, pp. 817-824. 

21. Pitz, R. W.; and Daily, J. W.: Experimental Study of Combustion 

in a Turbulent Free Shear Layer Formed at a Rearward Facing 
Step. AIAA Paper 81-0106, Jan. 1981. 


16 



Lewis motion-picture film supplement C-308 is available on loan. Requests will be filled in the 
order received. 

The film (16 mm, 15 min, color, sound) shows high-speed color schlieren motion pictures of a 
premixed propane-air flame stabilized by a rearward-facing step. A numerical technique that uses a 
random walk method for turbulence modeling and does not rely on empirical closure models was 
developed for analyzing turbulent combustion. Computer-produced color motion pictures of the 
volume fraction burned are used to compare the predictions with the experimental records. 

Requests for film supplement C-308 should to be addressed to 
NASA Lewis Research Center 
Attn; Chief, Management Services Division (60-1) 

21000 Brookpark Road 
Cleveland, OH 44135 
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